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Abstract. When a second-order phase transition is crossed at fine rate, the evolution 
of the system stops being adiabatic as a result of the critical slowing down in 
the neighborhood of the critical point. In systems with a topologically nontrivial 
vacuum manifold, disparate local choices of the ground state lead to the formation of 
topological defects. The universality class of the transition imprints a signature on the 
resulting density of topological defects: It obeys a power law in the quench rate, with 
an exponent dictated by a combination of the critical exponents of the transition. 
In inhomogeneous systems the situation is more complicated, as the spontaneous 
symmetry breaking competes with bias caused by the influence of the nearby regions 
that already chose the new vacuum. As a result, the choice of the broken symmetry 
vacuum may be inherited from the neighboring regions that have already entered 
the new phase. This competition between the inherited and spontaneous symmetry 
breaking enhances the role of causality, as the defect formation is restricted to a fraction 
of the system where the front velocity surpasses the relevant sound velocity and phase 
transition remains effectively homogeneous. As a consequence, the overall number of 
topological defects can be substantially suppressed. When the fraction of the system 
is small, the resulting total number of defects is still given by a power law related to 
the universality class of the transition, but exhibits a more pronounced dependence 
on the quench rate. This enhanced dependence complicates the analysis but may also 
facilitate experimental test of defect formation theories. 
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1. Introduction 

A phase transition is a transformation between different states of matter, and is typically 
induced by the quench of an external control parameter A through a critical point 
A c . The Kibble-Zurek mechanism (KZM) is a theory designed to describe the non- 
equilibrium dynamics in a scenario of spontaneous symmetry breaking and was originally 
developed for classical and continuous phase transitions [1, 2, 3, 4]. 

The original question of what happens in rapid phase transitions arose in cosmology, 
where it was realized that, as a result of relativistic causality, essentially all field theories 
predict formation of topological defects, with dramatic consequences for the subsequent 
evolution of the Universe [1, 2]. It was later proposed that a similar mechanism of 
topological defect formation should occur in all phase transitions traversed at the finite 
rate (including condensed matter phase transitions), although (instead of relativistic 
causality) the resulting density of defects should be tied to the critical scalings [3, 4]: 
The critical slowing down means that the newly forming phase will be able to coordinate 
the choice of the symmetry breaking over domains of limited size. We refer the reader 
to the earlier reviews [5, 6, 7] of the subject and recall the basics of the mechanism 
only briefly. When a system is quenched through a critical point of a second-order 
phase transition both the equilibrium correlation length £ and relaxation time r diverge 
according to 

f (*) = rhri r W = r^^> (!) 

in terms of the reduced parameter 

e=^. (2) 

Above, v and z are the critical exponents that define the universality class of the 
transition. The velocity with which the order-parameter phase information propagates 
(e.g. the speed of spin waves or second sound) is of the order of and no larger than the 
corresponding characteristic velocity 

s(t) = ^ = ^\e\t*- 1 ». (3) 

T[t) T 

Assume a linear quench e(t) = t/rQ, which drives the system through the critical 
point at t = 0. KZM simplifies the dynamics by distinguishing three stages: an adiabatic 
stage far away from the critical point, an impulse stage in the neighbourhood of the 
transition where the system is effectively frozen due to the divergence of the relaxation 
time, and a last stage governed by adiabatic dynamics again, sufficiently far away from 
the critical point and on the other (broken-symmetry) side of the transition. These 
three regimes are schematically represented in Fig. 1. Within this picture, the freeze- 
out occurs at the instant t when the relaxation time of the system r[e{t)] equals to 
the time scale of the quench, e/e. At that instant the dynamics approximately freezes 
(enters the impulse stage), and there is a breakdown of adiabaticity. Setting r(t) ~ |^/^|, 
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Figure 1. Impulse approximation for a second-order phase transition is central to the 
Kibble-Zurek mechanism. The equilibrium relaxation time r diverges at the critical 
point, here reached at t = via a linear quench e = t/rQ. When the time remaining 
for the transition becomes comparable to the relaxation time, the system cannot follow 
the quench and the dynamics is no longer adiabatic - it is no longer in equilibrium 
defined by the momentary value of e(t) (shaded area). The Kibble-Zurek mechanism 
estimates the typical size of the domains in the broken symmetry phase by the value 
of the equilibrium correlation length at the freezeout time t. 



it follows [3, 4] that 

£~ (tot^)^. (4) 

The key insight of the KZM is that the average size £ of the domains in the broken 
symmetry phase is set by the equilibrium correlation length evaluated at the freezeout 
time. This appears to hold even for transitions where only one side in Fig. 1 is effectively 
present [8, 9]. 

A universal power-law scaling of the size £ of domains that can coordinate 
spontaneous choice of broken symmetry on the quench rate emerges as a result, with 
the power given by a combination of the critical exponents, 

V 

i = at) = eo (^) ^ . (5) 

The symmetry in Fig. 1 suggests that freezing occurs at t ~ — £, and that any dynamics 
in the subsequent impulse stage, for t £ [— £, t], can be disregarded. However, numerical 
studies with asymmetric quenches around the critical point indicate that both sides may 
matter, and suggest that the domain formation may be decided at t after the transition, 
at t « +t. [10, 8, 11, 9]. The density of excitations resulting from the quench scales 
as d ~ where D is the dimensionality of the system. For topological defects this 
means a single "piece" of the defects per domain of size £. 
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The prediction of defect density based on Eq. (5) suggests that the experimental 
study of non-adiabatic dynamics across a second-order phase transition may be used 
to probe the universality class of a system. In particular, one might be able to use it 
to determine experimentally the dynamical exponent z which is often more difficult to 
measure than v. 

The adiabatic-impulse approximation can be invoked to describe the dynamics of 
a Landau-Zener transition: KZM scaling yields surprisingly accurate predictions [12], 
and can be related to avoided level crossing scenarios [ ] (e.g., so that the asymmetry 
between the two sides of the transition can be also probed) . The KZM is also applicable 
to quantum phase transitions (QPTs) - abrupt transitions between the nature of the 
ground state of quantum many-body systems that occur at zero temperature as a 
consequence of the gradual change of its Hamiltonian [14, 15, 16, 17]. The control 
parameter A can be now related to e.g. an external magnetic field or to the strength 
of the interactions. The transition is predicted to leave the system with the density 
of excitations given by KZM scalings (see reviews [ , ]). The relaxation time is then 
governed by the inverse of the energy gap between the ground state and the first excited 
state. This gap closes according to A <~ |A — A c |^, and the original reasoning that yields 
i and £ can be deduced. In addition to the resulting expected value of the defect density 
the amount of entanglement [ ] as well as the effectiveness of the many-body system 
as a decohering environment [ ] can be also deduced using KZM paradigm in at least 
some cases. 

2. The inhomogeneous Kibble-Zurek mechanism 

The Inhomogeneous Kibble-Zurek mechanism (IKZM) is an extension of the 
paradigmatic KZM that allows one to describe the dynamics of a phase transitions 
in systems under an inhomogeneous quench A(r,t), or with a spatially varying critical 
point A c (r), or a combination of both scenarios. 

The earliest study of the effect of an inhomogeneous quench [ ] was primarily 
in the context of defect formation in 3 He. Experiments had been performed on phase 
ordering produced by thermal neutron injection into superfluid 3 He-B [21, 22]. When 
a neutron is absorbed via the exothermic reaction n + 3 He — » p + 3 H + 0.76 MeV, 
a small volume is raised above the critical temperature and as it subsequently cools a 
network of vortex lines is formed. Clearly the transition here is inhomogeneous; there 
is a thermal gradient and therefore the second-order transition occurs as a propagating 
front. The key parameter is the velocity vp = { T Q\^ e \)~ l °f the front. When it is large 
the effect is very similar to the homogeneous case. But when it is slow the phase of the 
order-parameter is constrained by that in the ordered state behind the front, so defect 
formation is restricted. Because of critical slowing down, there is always a region near 
the front where s < Vp : where s is given by Eq. (3), and therefore defects will form. In 
the case of neutron-induced defect formation in 3 He-B, it turns out that this region is 
of the same size as the bubble heated to above the critical point, so the homogeneous 
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Figure 2. Formation of solitons during the BEC induced by a thermal quench 
of a inhomogeneous atomic cloud: a) Isodensity contour of the gas in the harmonic 
trap. The critical temperatures for Bose-Einstein condensation depends on the density, 
and is highest in the centre of the trap. That is where the condensation will start. 
In sufficiently slow quench the phase selected by that central region would spread 
throughout the trap, resulting in a an equilibrium pool of BEC. b) However, for 
faster quenches, the choice of the ground state in the broken symmetry phase differs 
in different parts of the system and topological defects can form. Here, in the 
approximately one-dimensional "cigar" these are (semi-) topological grey solitons that 
correspond to sharp changes of the phase of the condensate. They will form only near 
the center of the trap, where the speed of the phase front exceeds the speed of the 
sound, and the transition is in effect homogeneous. As the phase front slows down 
near the ends of the "cigar", speed of sound may exceed the speed of the front, and 
the phase of the last region that selected the phase spontaneously will be then passed 
on (much like the preferred lattice orientation in the crystals grown by Czochralski 
process). Reprinted with permission from [ ]. Copyright 2009 American Physical 
Society. 

KZM estimate of defect density should still be reasonable, as indeed was confirmed by 
the experiments. On the other hand, in 3 He-A the transition is much slower and so 
fewer defects are formed. Moreover, as we show below, the dependence of defect density 
on quench rate is much stronger in this case. 

The approach below was developed for the study of soliton formation in the course 
of Bose-Einstein condensation (BEC) [ ]. Gaseous BECs are confined in traps that 
result typically in a "cigar" shape (see Fig. 2) of the resulting condensate, which is 
approximately one-dimensional, and is densest in the center of the cigar (where BEC 
forms first). 

Let us consider the case when A c = A c (x) (for simplicity we assume a one- 
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dimensional system) and the external quench is homogeneous and linear in time, 



\{t) 



Ao ( 1 - - 



(6) 



for t G [—tq, tq]. The choice A = A c (0) will be useful in the following (any other choice 
can be cast in this form by redefining t'q = tqA c (0)/Ao). 
The reduced control parameter becomes 
X c (x) - X(t) 



e(x, t) 



X c (x) 



This leads to an effective quench rate with a spatial dependence, 



t q( x ) = 



de(x, t) 



dt 



(7) 



(8) 



One can identify the location of the front xp crossing the transition at tp from the 
condition e(xp,tp) = 0. In particular, criticality is reached at point xp = x at a time, 

A c (x 



tp(x) = Tq 



1 



A c (0) 



(9) 



(10) 



This allows us to conveniently rewrite Eq. (7) as 
e(llt ) = *£Wzi, 

which resembles the expression in the homogeneous scenario. We may pursue three 
alternative routes to identify the characteristic value of t) which determines the size 
of the domains in the broken symmetry phase. 

• One way to go is to invoke the adiabatic-impulse-adiabatic approximation to 
determine the value e(x,t) at which the system departs from the instantaneous 
equilibrium. It suffices to match 



r[e(x,t)] = = £ ( x > t ) T Q( x )> 

t) 



where the dot means time-derivative with respect to T = tp — t. It follows that 



e(x) 



T"0 



l + vz 



(12) 



• As an alternative, one can identify six) by looking at the rate of change of the local 
equilibrium correlation length 



v 6 (e) 



d£, de 
dedi 



(13) 



T Q {x)\e{x,t)\^- 

required to match the velocity of perturbations, the relevant sound velocity which 
can be upper bounded by the ratio of the local correlation length and the relaxation 
time, 



\u(z-l) 



(14) 



Solving s[s(x,t)] = V£[e(x,t)] one finds e(x) = (y 



tq(x)< 



l + vz 
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A third alternative is based on the comparison of the local spatial extension of the 
sonic horizon h{x) with the equilibrium correlation length. Here 



h[e(x, t)] 



s{T f )dT f 



s[e{x ) t)}e{x ) t)T Q {x) 



D + i 



(15) 



[u(z - 1) + 1] 



t q( x ) 



Solving h[e{x,t)\ = £[s(x,t)] we find £{x) 

Up to the prefactors given by the critical exponents, these alternative approaches 
agree on the prediction for e(x), associated with the local freeze-out time T{x) = 
[ t o t q( x Y z ] 1+uz • It follows that the typical size of the domains in the broken symmetry 
phase is given by 

~ t q( x )~ 



TO 



l + vz 



(16) 



This prediction is consistent with the use of a local density approximation for tq(x) in the 
KZM length scale predicted for homogeneous systems. However, the role of causality is 
enhanced during the defect formation across a inhomogeneous phase transition. Defects 
are expected to be formed only where the velocity of the front surpasses the characteristic 
velocity of perturbations [20, 24]: in that case, the transition is effectively homogeneous, 
in that different fragments of the condensate cannot pass on their local choice of 
broken symmetry to their neighbours. Numerical and analytic studies in specific models 
[24, 25, 26, 11, 27] have confirmed this expectation (see e.g. Fig. 3). 
The speed of propagation of this front can be estimated to be 



v F 



dxi 



dti 



A c (0) 



d\ c (x) 



dxi 



(17) 



As for the characteristic local velocity of the perturbations, it is given by the second- 
sound velocity that can be upper bounded by the ratio of the local frozen out correlation 
length and the relaxation time scale r{x) — r[s] = this is, by 



TO 



tq{x) 



1+vz 



(18) 



The following cases can be distinguished. If v F < s everywhere in the system, the 
resulting dynamics is adiabatic and defect formation is suppressed. As a result, IKZM 
may offer a way to an adiabatic dynamics across a second order phase transition. 

The sharpness of this inequality has been investigated in a variety of models 
including the stochastic Ginzburg-Landau equation [24], as seen in Fig. 3a. It has also 
been studied in the ID quantum Ising (see Fig. 3b) and XY models [25, 26], the Gross- 
Pitaevskii equation for the miscible-immiscible phase transition of a two-component 
BEC [ ], and the Langevin dynamics of underdamped charged particles [ ]. When 
Vp > s everywhere in the system, the dependence of the density of excitations on the 
quench rate is correctly described by the homogeneous KZM. Clearly, this is always the 
case in homogeneous systems where vp diverges given that X c (x) = A c (0). 

A novel scaling results when vp > s holds only on an effective fraction X of the 
system, depending explicitly on the quench rate. Depending on the topology of the 
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Figure 3. Number of defects as a function of the front velocity: a) The 
simulations correspond to the dynamics of a one-dimensional real scalar field <j)(x, t) 
whose dynamics is governed by a time-dependent Ginzburg-Landau equation of the 
form <9 2 (/> + rjd t (j) — <9 2 (/> + [0 3 — e(x,t)4>] = 0(x,t). Here, 6(x,t) is a Gaussian white 
noise satisfying (6(x,t)6(x' ',£')) = 2r]6oS(t — t')5(x — x') and the inhomogeneous 
quench describes a sharp pressure front e(x,i) = eosgn(t — x/vp) propagating at 
velocity vp = v. The role of s is played by the characteristic sound velocity s = 
[1 + 2?7/3y / eo] -1 / 2 , which takes the value s = 0.83 for r\ = eo = 1. From top to bottom 
the curves correspond to the temperature #0 = 10 _1 , 10 -2 , 10 -4 , 10 -6 , 10 -8 , 10 -10 
and show that at low temperatures the onset of excitations is well captured by the 
inequality vp > s [ ]. Copyright 1999 by the American Physical Society, b) 
Formation of excitations in the inhomogeneous ID quantum Ising Model described 
by the Hamiltonian H— = — Yln=i Qn&n ~ ^2n=i a n a n+i- ^ n ^ ne homogeneous case 
(9n — g) m the thermodynamic limit, the system has critical points at g = ±1 
separating a paramagnetic phase (\g\ > 1) from a ferromagnetic phase (\g\ < 1). 
The simulation corresponds to a chain of N = 400 spins undergoing an inhomogeneous 
quench of the magnetic field of the form g n = 1 + tanh[a(n — vt)], driving locally the 
chain from the paramagnetic (g = 2) to the ferromagnetic phase (g = 0). In this 
system s = 2. The formation of quasiparticle excitations is suppressed for v < s. 
Adapted with permission from [25]. 



system there might be several disconnected regions where this occurs [28] . The number 
of excitations is then approximately given by 

rf— y / dx-^- T ° 1+VZ . (19) 

L J{x\v F >s} £0 Lrg(x)J 

The condition vp > s defines the effective size of the system where defect formation 
is allowed. Let us assume this is a simply connected region [—£,£]. Generally, the 
expression in Eq. (19) does not lead to a power-law for d in tq, and exhibits a more 
complex behavior dependence [28]. However, when the fraction X = x/L <C 1 it is 
possible to derive analytical scaling power law. In several instances (including the case 



Causality and non- equilibrium second- order phase transitions in inhomogeneous systems^ 



of Bose-Einstein condensation illustrated in Fig. 2) it has been found that [23, 29, 30, 28] 
d^J^ocf^-)^. (20) 



£(0)L \tq, 

For mean-field critical exponents {y = 1/2, z = 2) this results in a power-law with 
= 1, four times larger than in the KZM where the exponent is = 1/4. This 
pronounced dependence on the quench rate is particularly amenable to tests in the 
laboratory. 

The IKZM applies as well to quantum phase transitions [31, 25, 26, 6, 7]. The most 
important effect is similar to what we have seen in the classical phase transitions: Once 
the velocity of the front falls sufficiently below the local relevant sound speed, defect 
production is suppressed. Dziarmaga and Rams [25, 26] have provided an especially 
convincing discussion of this phenomenon in the quantum Ising and XY models in one 
dimension. 



3. Laboratory systems 

3.1. Trapped ions 

Experimental evidence supporting the IKZM was recently obtained in a structural phase 
transition in ion Coulomb crystals. Trapped ions constitute a versatile platform for 
simulating many-body physics [32]. The formation of structural defects in ion crystals 
by quenching the external potential was proposed in [29, 30]. An ion chain confined in 
a highly anisotropic Paul trap undergoes a continuous phase transition from a linear to 
a (doubly degenerate) zigzag structure when the transverse trap frequency is quenched 
across the critical value, see Fig. 4a. The interaction potential is given by 

1 N 1 2 i 

j = l j J 

for TV ions of mass m and charge e. The order parameter corresponds to the transverse 
size of the zigzag structure and obeys a Ginzburg-Landau potential [ ] . Laser cooling 
of the ions provides the relaxation mechanism by which vibrational excitations are 
damped. The possibility of testing the IKZM in this system arises from the axial 
harmonic confinement of frequency v ) which induces a modulation of the linear ion 
density of the form n(x) = n(0)(l — x 2 /L 2 ), where L is half the length of the chain and 
n (0) = tr* -^ or homogeneous system, the critical value of the transverse confinement 
is given by v±_ ~ 2. 051^0 where the frequency characterizing the Coulomb coupling 
reads ujq = [ 4^°^ ]"^ 2 - Using a local density approximation, this critical frequency 
acquires a spatial modulation v±_{x) = v±_[n(x)], making the phase transition ultimately 
inhomogeneous, with the quadratic inhomogeneity analogous to the one illustrated in 
Fig. 2 (and, as in Fig. 2, resulting from a harmonic trapping potential). A non-adiabatic 
crossing results in the formation of structural defects, kinks and anti-kinks, as shown 
in Fig. 4b. Numerical simulations in [29, 30] show that this non-equilibrium dynamics 
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(c) 




(b) 




Figure 4. (a) Structural phase transition of an ion Coulomb crystal in a harmonic 
trap. The chain consists of 30 172 Yb + ions in a linear Paul trap. As the transverse 
confinement is weakened there is a transition from the linear configuration to a 
doubly- degenerated zigzag chain, (b) Snapshots of an ion chain supporting a kink 
and an anti-kink, identified in the shaded areas, (c) Sequence of snapshots showing 
the inhomogeneous nature of the transition. As a result of the longitudinal axial 
confinement the inter-ion spacing is inhomogeneous. The transition is first crossed 
in the center of the chain, where the inter-ion Coulomb repulsion is higher. A front 
reaching the critical point spreads then sideways along the chain leading to the growth 
of the broken-symmetry zigzag phase. Courtesy of T. E. Mehlstaubler. 



is well described by the KZM. As the quench rate decreases, there is a crossover where 
among different power-laws, see Fig. 5. For sufficiently fast quenches the whole ion 
chain is effectively homogeneous and the effect of the harmonic trap can be ignored. As 
a result the typical size of the zigzag domains scales in the quench rate in agreement 
with the KZM, i.e. Eq. (5). As the quench rate decreases, it is possible to probe the 
inhomogeneous features of the system, and the scaling in Eq. (20) governs the defect 
density. The recent experiment at Physikalisch-Technische Bundesanstalt (PTB) [27] 
reported the scaling of defects as a function of the quench rate in an underdamped ion 
chain, which admits a description in terms of a Ginzburg-Landau theory with critical 
exponents v — 1/2 and z — 1. Typical realizations led to {0, 1} defects, as shown in 
Fig. 4. In that case of either or 1 defects (of arbitrary topological charge) in the 
whole system, in a homogeneous phase transitions a doubling of the KZM scaling was 
experimentally reported in [34]. Subsequent theoretical works predicted a doubling of 
the scaling in a variety of systems [35, 8, 36]. The PTB experiment found that a similar 
doubling occurs as well in the IKZM, and verified that the probability of forming a 
single defect pi scales as 

2(l + 2zy) 



Pi 



2x 1 2 I Tq\ 1+1/2 / T 



^fj - U 1 ■ (22) 



-«(0). 

where the last equality holds for v — 1/2 and z — 1, see Fig. 5. This result applies 
whenever 2x ~ £(0). Considering that the order parameter takes a random walk of 
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Figure 5. Scaling of the number of kinks in the quench rate for the linear to zigzag 
phase transition in an under damped harmonically trapped ion chain. The different 
lines are guidelines to the eye illustrating the crossover among the different power laws 
observed. At fast quench rates (right part of the plot) the system behaves as affectively 
homogeneous and the standard KZM prediction applies. The power-law is described 
by the inverse of £ in Eq. (5) with an exponent 1/3 for v = 1/2, z = 1. At intermediate 
quench rates causality determines the effective size of the system leading to a more 
pronounced dependence on the quench rate, as predicted in [29, 30]. In this regime the 
inhomogeneous Kibble-Zurek dictates a power-law exponent 4/3. At slower quench 
rates, the KZM correlation length £ becomes comparable to the effective system size 
~ 2£, and typical realizations leads to one or no defects at all. The density of defects 
exhibits then a doubling of the IKZM scaling [27] with a power-law exponent 8/3. This 
is the regime which would be probed experimentally in [27]. The mean kink number 
n was obtained by averaging over 2000 realizations the number of kinks per cycle in a 
N = 30 ion chain. Adapted from [27]. 

n 2x/£(0) steps, we expect a change of the behavior from the power-law in Eq. (22) to 
that of Eq. (20) as the number of steps increases from 0{1). Indeed, Langevin dynamics 
simulations show that in larger ion chains undergoing faster quenches multiple defects 
are formed, paving the way to the observation of the IKZM scaling in the absence of 
doubling. 

3.2. Bose-Einstein condensates 

Bose-Einstein condensation of trapped ultracold gases constitutes another relevant 
scenario where the inhomogeneous density of the cloud induced by the confinement 
is expected to affect the non-adiabatic dynamics. Much intuition about the transition 
can be gained by considering the Gross-Pitaevskii energy functional describing a trapped 
(scalar) Bose gas, which is given by 

£ GP [$] = j (|^l w l 2 + [ y W - #l 2 + \ l $ l 4 ) rf3r > ( 23 ) 
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Figure 6. Spontaneous formation of vortices during the growth of a Bose-Einstein 
condensate, a) To illustrate the growth dynamics of a BEC in a 3D harmonic trap, 
snapshots of three different expanded BECs are obtained at different points in the BEC 
growth cycle. Note that these images do not show the growth of a single BEC. In left 
picture a few small regions of high density appear to stand out against a thermal-cloud 
background, although it is possible that these only appear as separated regions due to 
imaging artifacts. The central picture shows a BEC in an early stage of growth, with 
a single clear vortex core and a few smaller dark regions that may represent additional 
vorticity. In the snapshot at the right, a single vortex is clearly seen in a fully formed 
BEC. All images taken after expanding from the harmonic trap, b) Merging of three 
independent Bose-Einstein condensates might lead to the spontaneous formation of 
a vortex, whenever the circulation of the phase around a given point adds up to a 
multiple of 2tt. Courtesy of B. P. Anderson. 

where \i is the chemical potential and g in the interaction coupling. For a homogeneous 
system the external potential can be set V(r) = 0, and Bose-Einstein condensation 
is expected whenever \i > 0. This transition corresponds to a U(l) spontaneous 
symmetry breaking. In trapped systems, one can introduce the effective chemical 
potential /x(r) = \i — V(r). Nonetheless, the measurements of the critical exponent 
v indicate deviations from mean-field theory [ ]. At the U. of Arizona, experiments on 
forced evaporative cooling of a thermal cloud through Bose-Einstein condensation have 
demonstrated the spontaneous formation of vortices in harmonic and toroidal traps [38] . 
As the condensate starts to grow, the phase of the condensate wavefunction is chosen 
in patches of area ^ £( r ) 2 • As these patches merge with each other, there is a chance 
for the circulation of the phase around a given point to become a multiple of 2n and 
lead to the formation of a vortex. This is illustrated in Fig. 6a where three different 
stages of the growth of a BEC are shown. Indeed, the Tucson group demostrated in 
a different experiment that the explicit merging of three independent BECs does lead 
to the formation of vortices [ ] (see Fig. 6b), with a probability determined by the 
geodesic rule [ ]. In [38], the limitation of the achievable quench rates and available 
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statistics prevented the test of the KZM scaling. It was later pointed out that the 
experimental setup could explore both homogeneous and inhomogeneous scenarios [28] . 

When the atomic cloud undergoing BEC is highly anisotropic, soliton creation 
follows from the non-adiabatic dynamics, and the dependence of the number of 
spontaneously formed solitons on the quench rate observed in the homogeneous scenario 
[41] can be dramatically affected by the presence of the trapping potential [23]. 
Numerical simulations using the classical field method exhibit a power-law scaling [42] in 
agreement with the IKZM scaling predicted in [ ] . To assess the observability of IKZM 
during Bose-Einstein condensation, the above studies should be supplemented with a 
kinetic theory of evaporative cooling, relating the quench rate of the forced-evaporation 
to the effective cooling rate of the cloud. 

Inhomogeneities might be of relevance to the dynamics of other ultracold atom 
systems, as in the phase separation and pattern creation in a binary Bose-Einstein 
condensate [43, 11], the spontaneous spin vortex formation in a trapped ferromagnetic 
Bose-Einstein condensate [4 ], or the quantum phase transition from Mott insulator 
to superfluid phase, where a power-law scaling for the density of excitations has been 
experimentally demonstrated [45]. A universal scaling of the density of excitations 
governs as well the coherent quantum evolution of a many-particle system following a 
slow amplitude quench of a confining potential [46]. 

In experiments with ultracold gases and trapped ions the spatial dependence of 
A c (x) often arises as a result of the density modulation induced by an external trapping 
potential. We close this section noticing that it is possible to introduce a trap-size 
scaling of the correlation length with the trap size [ ] . This finite-size scaling modifies 
the divergence of the correlation length in the neighbourhood of the critical point. The 
discussion of the IKZM above is based on the power-law divergence of the equilibrium 
correlation length for a homogeneous system and is therefore restricted to quenches 
where the value of £ is small compared to the trap size. 

4. Defect formation beyond scaling laws 

In this section we look for a more general expression for the formation of kinks, removing 
the restriction I < 1. We shall derive a general expression for the density of defects 
which cannot be cast as a power-law of the quench rate. Nonetheless, we shall recover 
the results of the previous section for X <C 1. 

For the sake of illustration, we shall focus on systems where the spatial dependence 
of the critical control parameter is described by a Thomas-Fermi profile, i.e. 

A c (x) = A c (0)(l-X 2 r, (24) 

where a > 1 and X = x/L. Particular cases correspond to the linear-to-zigzag transition 
in harmonically trapped ion chains [29, 30], the mixing- demixing transitions [43, 11], or 
the magnetic transition in spinor BEC [ ]. A Gaussian spatial dependence of A c (x) is 
of relevance to non-interacting ultracold gases and was discussed in [23, 28]. Consider 
a linear quench of the form A(i) = A C (0)(1 — £/tq), for t G [— tq, tq]. 
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We note that \ c (x) reaches its maximum at x = so the critical point is first 
reached at the center. Symmetry breaking spreads sideways, and the front crossing the 
transition reaches point a; at a time 

X c (x)' 



t, 



tq[1 



A c (0) 



Relative to this time, it is possible to rewrite 

where the local quench rate is 

X c (x) 



t q( x ) = T Q 



A c (0) 



r Q (i-xr- 



(25) 



(26) 



(27) 



For the Thomas- Fermi profile the KZM bound to the velocity of perturbations reads 



r \ T Q / 



while the front velocity is given by 
L 1 



(28) 



(29) 



T Q 2a\X\(l -X 2 )"" 1 ' 

Near the origin X w 0, v-p diverges and defect formation is expected. Indeed, the 
inequality wf(^f) > s allows us to find that defect formation is possible for all 
, where 

L 2 



x 



TO 



l + vz 



2a£ \tq 
The number of defects is given by 



0(X 6 ). 



(30) 



(I ~ y 



dx — 

-x Co 



TO 



l + vz 



2x 



tq(x) 
ctv 13 
T+vz' 2' 2' 



(31) 



where 2 ^i(^,6, c\z) is the hypergeometric function, provided that x < L, as is the 
case. As an upshot, whenever the spatial dependence of the effective quench rate tq(x) 
becomes important over the domain where defect formation is allowed by causality, 
the resulting number of kinks is not described by a power-law scaling any more. Slow 
quenches and tight traps provide actually the ideal scenario to observe IKZM scaling, 
see [ ], end of Sec. 3 (assuming no decay mechanisms). For x«L, 



d 



2x 



1 - 



OLV 



X 



3(1 + vz) L 2 
0(X 3 ) 



+ 0(X A ) 



(32) 



l+2v 
l + vz 



<o \tq, 

where we recognize the leading term as the estimate for the density of defects in the 
previous subsection. Whenever x can be approximated by Eq. (30), this term leads to 



(33) 
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a power-law scaling as a function of the quench rate. One can identify two sources of 
deviations from this scaling. First, the exact solution of ^f(^f) > s for x, which typically 
leads to a dependence on tq which is not a power-law. As mentioned above, it can be 
the case that ^f(^f) > s holds in spatially separated regions of the system leading to 
a more complicated dependence of X on tq [28]. Second, the local dependence of the 
quench rate tq(x), which is responsible for the higher leading terms in the expansion in 
brackets, Eq. (32). 

So far, we have restricted the presentation to one-dimensional systems. However, 
the extension to higher dimensions is straightforward. For the sake of illustration, 
consider a two dimensional systems with radial symmetry, area A = 7rL 2 , and critical 
control parameter A c (r) = A C (0)(1 — R 2 ) a where R = r/L. The density of defects of size 
£(t) 2 is given by 

1 f r 2nrdr ( r \ 1+vz 

2v 



kL 2 J Co \tqM, 
1 fr \^ \-{l-R 2 f 



^ VqJ p 

= d KZM f(R). (34) 

Here, /3 — 1 — and we have assumed that the causality condition restricts the 
formation of defects to r G [0,f) and defined R — r/L. This expression is particularly 

interesting since the term d^zM = -&[ — ) can be recognized as the KZM estimate 
for the density of defects in a homogeneous system. The role of the inhomogeneity 
is encoded in the additional term f(R), which can be expanded whenever R <C 1 to 
f(R) = R 2 + 0(R 4 ). Noting that f can be approximated by the form of x in Eq. (30), 
we derive the power-law scaling 

2(l + 2i/) 

rf = A " +0{R A ), (35) 

2q £o \ t qJ 

For example, this dependence could be experimentally probed by studying vortex 
formation during Bose-Einstein condensation induced by a slow thermal quench in a 
tight trap [28]. 

We note that in these higher dimensional systems there exist the possibility for 
the critical control parameter to be inhomogeneous only along a given set of degrees 
of freedom. This results in power-law scalings with exponents intermediate between 
those of the homogeneous and inhomogeneous scenarios. An example of relevance to 
the experiment [3£ ] is the spontaneous vortex formation during the growth of a Bose- 
Einstein condensation in toroidal trap [28] . 

5. Discussion: Inhomogeneous driving as a tool to suppress excitations 



Suppression of excitations across a phase transition is of interest in a wide variety 
of applications, including the preparation of quantum phases in the ground state 
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in a quantum simulator [48, 49], and implementing protocols for adiabatic quantum 
computation [50]. Different approaches have been put forward to approach or mimic 
adiabaticity. Finite size-effects pave the way to an adiabatic dynamics, whenever the £ 
is larger than the system size L. An adiabatic dynamics can be implemented exploiting 
the energy gap that arises from the finite size of the system [ ] . An alternative approach 
is to consider a finite-rate non-linear quench of the form A(i) = A |t/rg| r and optimize 
the power-law exponent r to minimize the density of excitations following the quench 
through the critical point [52, 53]. The optimal power turns out to be given by the 
logarithm of the quench time multiplied by universal critical exponents characterizing 
the phase transition [ ]. However, an accurate knowledge of the critical point and 
excellent control of A(i) are of essence in this approach. Otherwise, the quench might be 
linearized in the proximity of the critical point. The KZM scaling can be modified as well 
by coupling a quantum critical system to an environment. The resulting nonequilibrium 
dynamics exhibits corrections to the KZM scaling and paves the way to relaxation of 
excitations [ ]. The topology of the system can influence production of defects. This 
fact was illustrated in the Creutz model, which supports topological edge states [55]. 
It was shown that the dynamics following a quench of an edge states dramatically 
increases the exponent governing the power-law scaling and diminishes the density of 
excitations. This observation was extended in a subsequent study of the adiabatic 
dynamics of Majorana fermions across a quantum phase transition. It was shown 
that the KZM describes the formation of bulk defects but it is not valid for for edge 
Majorana fermions. The root of these non-universal deviations was pinned down in the 
localization dynamics of edge states [ ] . Further advances in tailoring excitations have 
been achieved using optimal quantum control (OQC) strategies. The combination of 
time-dependent Density Matrix Renormalization Group (t-DMRG) with the Chopped 
RAndom Basis (CRAB) optimization algorithm allows to tackle a great variety of non- 
integrable quantum many-body one-dimensional systems [57, 58]. Remarkably, OQC 
allows one to optimize adiabatic evolution with neither knowledge of the instantaneous 
eigenstates nor the resulting gap of an arbitrary time-dependent system [ ]. Of course, 
the applicability of OC in many-body systems is not restricted QPT, and can be applied 
to design quenches without crossing a critical point [ ]. Yet another approach, consists 
in using counterdiabatic interactions which allow to drive the dynamics through the 
instantaneous eigenstates of the many-body system [61]. 

The IKZM can be though of as a complementary tool to these various approaches 
since it helps to reduce the number of excitations generated during the crossing of 
a phase transition [ ]. Its relevance to many-body state preparation and adiabatic 
quantum computation was pointed out in [25, 26], see as well [ ]. This line of thinking 
has recently been applied to the adiabatic growth of an antiferromagnetic phase in a 
Mott insulator made of two atomic species [62] . 



Causality and non- equilibrium second- order phase transitions in inhomogeneous systems!! 
6. Conclusion 

The divergence of the relaxation time in the neighbourhood of the critical point induces 
a non-adiabatic dynamics in which causally disconnected regions perform independent 
choices of the broken-symmetry vacua. As result, topological defects are formed, their 
density scaling with the quench rate as described by the Kibble-Zurek mechanism. In 
inhomogeneous systems, the local choice of the ground state in the broken symmetry 
phase can be communicated to neighbouring regions. Criticality is reached locally as 
the phase transition front spreads gradually across the system. When its velocity is 
lower than the speed of the relevant sound the formation of defects is suppressed. 
Otherwise causality determines the effective size of the system. This competition 
between local, spontaneous choice of broken symmetry and the bias of the neighbouring 
regions that had already made their selection leads to a more pronounced dependence 
of the defect density on the quench rate than in the homogeneous scenario. The latter 
is only recovered in inhomogeneous systems for fast quench rates. Regimes where the 
dependence is captured by a simple power-law have been identified, as well as a variety 
of crossovers among different power-laws. This universal nonequilibrium dynamics is of 
relevance to a wide variety of laboratory systems including trapped ions and ultracold 
gases. 
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